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Abstract. The relationships between the effective surface 
(Teff) and internal temperatures of neutron stars (NSs) 
with and without accreted envelopes are calculated for 
Teff > 5 X 10^ K using new data on the equation of state 
and opacities in the outer NS layers. We examine various 
models of accreted layers (H, He, C, O shells produced by 
nuclear transformations in accreted matter). We employ 
new Opacity Library (OPAL) radiative opacities for H, 
He, and Fe. In the outermost NS layers, we implement the 
modern OPAL equation of state for Fe, and the Saumon- 
Chabrier equation of state for H and He. The updated 
thermal conductivities of degenerate electrons include the 
Debye-Waller factor for the electron-phonon scattering in 
solidified matter, while in liquid matter they include the 
contributions from electron-ion collisions (evaluated with 
non-Born corrections and with the ion structure factors 
in responsive electron background) and from the electron- 
electron collisions. For Tea < 10®'^ K, the electron con¬ 
duction in non-degenerate layers of the envelope becomes 
important, reducing noticeably the temperature gradient. 
The accreted matter further decreases this gradient at 
Teff > 10® K. Even a small amount of accreted matter 
(with mass ^ IO^^^Mq) affects appreciably the NS cool¬ 
ing, leading to higher at the neutrino cooling stage 
and to lower Teff at the subsequent photon stage. 

Key words: stars: neutron - pulsars: general - dense 
matter 


1. Introduction 

It is well known that the cooling of a neutron star (NS) is 
strongly affected by the relationship between the internal 
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temperature of the star, Tb, and its effective surface tem¬ 
perature Teff. Throughout this paper, Tb denotes the tem¬ 
perature at the outer boundary of the isothermal internal 
region. The Tb-Teff relationship is determined by equa¬ 
tion of state (EOS) and thermal conductivity of matter in 
the outer NS envelope. For non-magnetized NSs with en¬ 
velopes composed of iron, this relationship has been thor¬ 
oughly studied in a classical article of Gudmundsson et al. 
(1983, hereafter GPE). Several papers (Hernquist 1984, 
1985, Van Riper 1988, 1991, Schaaf 1988, 1990) have con¬ 
sidered strongly magnetized NS envelopes. 

In this article, we will reconsider thermal transport in 
non-magnetized envelopes of NSs and extend the results 
of the preceding authors in two respects. 

First, we analyze matter composed not only of iron, 
but of light elements as well. The light elements can 
be provided by an accretion from a supernova remnant 
(e.g., Chevalier 1989, 1996, Brown & Weingartner 1994), 
from interstellar medium (e.g., Miralda-Escude et al. 1990, 
Blaes et al. 1992, Nelson et al. 1993, Morley 1996), from a 
distant binary component, or by comets. Freshly accreted 
matter burns then into heavier elements (He, C, O, Fe) 
while sinking within the NS. Recent multiwavelength ob¬ 
servations of the Geminga pulsar (IE 0630-1-17.8) suggest 
a possible H or He cyclotron feature in its spectrum (Big- 
nami et al. 1996). If confirmed, it may be a direct obser¬ 
vational evidence of the presence of the light elements in 
the pulsar atmosphere. Chemical composition affects the 
EOS and thermal conduction, and, therefore, the thermal 
structure and cooling of a NS. As a reference case, we 
reconsider the outer NS envelopes composed of iron. 

Second, we implement new, advanced theoretical data 
on EOS and thermal conductivity of dense matter. Specif¬ 
ically, we employ the Opacity Library (OPAL) radiative 
opacities for H, He, and C, improved considerably with 
respect to the Los-Alamos opacities used in the previous 
studies. In the outermost NS layers, we also implement the 
modern OPAL EOS for Fe and the Saumon-Chabrier EOS 
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for H and He. We use improved thermal conductivities 
of degenerate electrons. For solidified matter, we employ 
the thermal conductivity due to the electron-phonon scat¬ 
tering obtained with the inclusion of the Debye-Waller 
factor. For liquid matter, we recalculate and implement 
the thermal conductivity due to Coulomb electron-ion and 
electron-electron collisions. The electron-ion scattering is 
described with the exact (non-Born) Coulomb cross sec¬ 
tions and with the ion structure factors calculated when 
taking into account the response of the electron back¬ 
ground. The new physics input allows us to extend the 
results of previous studies to colder NSs, with Tes down 
to 50000 K. 

The physics input is described in Sec. 2. In Sec. 3, we 
calculate the Tb-Jeff relationships for non-accreted and 
partly accreted NSs and analyze the sensitivity of the re¬ 
sults to the uncertainty in our knowledge of the electron 
thermal conductivity. In addition, we simulate the cool¬ 
ing of NSs with standard and enhanced neutrino energy 
losses. We show that the cooling of a NS with the ac¬ 
creted envelope can be quite different from the cooling of 
a non-accreted NS. This can change the conclusions on 
the internal structure of NSs deduced from comparison of 
theoretical cooling curves with observations of NS thermal 
radiation. Useful analytical formulae for the physics input 
are given in the Appendix. 


2. Physical input 

2.1. Thermal structure equation 

Consider thermal transport throughout the outer envelope 
of a NS that extends from the surface to the layers with 
the density pb ^ 10^° g cm“^. This envelope, of which 
study is the aim of the present paper, provides the main 
thermal insulation of the NS interior. The envelope is thin 
(~ 10^ m deep) and contains very small fraction (~ 10“^) 
of the NS mass (GPE). In these layers, one can neglect the 
nonuniformity of the energy flux due to the neutrino emis¬ 
sion and the variation of the gravitational acceleration. 
Then the temperature profile obeys the thermal structure 
equation (GPE), which can be written as (e.g., Van Riper 
1988) 


composed of the radiative opacity A'rad and the equivalent 
electron conduction opacity Kc. The latter is related to the 
electron thermal conductivity k as 


K, 


3pK ’ 


( 3 ) 


where a is the Stefan-Boltzmann constant. 

The effective temperature Teg is defined through the 
Stefan’s law. 


.^rad — 47ri?^(TT^, 


( 4 ) 


where £rad is the local radiative luminosity. The apparent 
luminosity measured by a distant observer is Coc = (1 — 
rg/R) £rad, and the apparent surface temperature inferred 
by the observer from the spectrum is Too = TeSy/l — r^jR 
(e.g., Thorne 1977). 

We adopt the standard (GPE) outer boundary condi¬ 
tion to Eq. (P by equating Tes to the temperature Tg at 
the stellar surface in the Eddington approximation. The 
radiative surface, determined in this way, lies at the opti¬ 
cal depth r = 2/3, that is approximately at 

Pg«(2/3)5/ifg. (5) 


Here and hereafter the subscript “s” denotes quantities 
taken at the surface, while the subscript “b” is assigned 
to those at the inner boundary p = pb- Previous calcu¬ 
lations of the radiative transfer in the outer part of the 
atmosphere (Zavlin et al. 1996) have shown that, for the 
parameters of interest, there is no convective instability at 
T < 1 which could invalidate the Eddington approxima¬ 
tion. 

Starting from the surface specified by Eq. (^, we inte¬ 
grate Eq. (H) inward, using the classical 4-step 4th-order 
Runge-Kutta algorithm (e.g., Fletcher 1988). The inte¬ 
gration step varies with depth and temperature to en¬ 
sure the variations of T and K to be smaller than 5% 
at one step. Following previous studies (e.g., GPE, Hern- 
quist 1985, Van Riper 1988), we terminate the integration 
at pb = 10^° g cm“^. A comparison with analytically solv¬ 
able models (e.g., Hernquist & Applegate 1984) shows that 
our numerical procedure determines Tb at a given Tg with 
an error ^1%. 


dlogr 3 PK T^^ 
dlogP 16 p T4 ■ 

Here, P is the pressure, Tgs is the effective temperature, 
g = GM/{R y/1 — rg/R) is the surface gravity, M is the 
stellar mass, R the stellar radius, G the gravitational con¬ 
stant, and Cg = 2GM/c^ is the gravitational radius of the 
star. Furthermore, K is the opacity. 


K = 


1 

R-rad 



( 2 ) 


2.2. Equation of state 
2.2.1. Plasma parameters 

Although the density p does not enter Eq. ( 0 ) explicitly, 
its value is required to determine the opacity K at given 
(T, P). Thus, we need an EOS of matter in the outer NS 
envelope. In this paper, we are specifically interested in 
H, He, C, O, and Fe plasmas, which are likely to compose 
partly accreted NS envelopes. The parameters of interest 
cover the region which extends from the partially ionized 
non-degenerate non-relativistic plasma near the surface. 
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through the domain of pressure ionization, to high-density 
layers of fully ionized relativistic plasma. We assume a 
strongly stratified envelope, i.e., that at each given den¬ 
sity, the plasma is composed of electrons and of ions of 
one chemical element which can be in different ionization 
stages. 

The plasma density can be characterized by the pa¬ 
rameter Tg = Oe/aB, where Oe = (47rne/3)“^/^ is the 
mean inter-electron distance, Ue is the number density 
of electrons, and ob is the Bohr radius. The parameter 
Tg is directly related to the relativistic parameter of de¬ 
generate electrons x = pF/rUeC = 1.009(p6(^)/(^))^^^ 
(e.g., Yakovlev & Shalybkov 1989): = 0.0140/a;. Here, 

Pf = is the electron Fermi momentum, me is 

the electron mass, pe = p/10® g cm“^, and {Z) and (H) are 
the mean charge and mass number of ions, respectively. In 
our case, the averaging is over different ionization stages 
of ions of the same chemical element. 

The electron degeneracy temperature is given by 


Tf = 



1 meC^/fcB, 


( 6 ) 


where fcs is the Boltzmann constant. 

Let us also introduce the ion coupling parameter, 
r = /{tteksT), which measures typical energy of 

ion Coulomb interaction relative to the thermal energy (e 
being the electron charge). 


2.2.2. High density regime 

Let us specify the high density regime as Tg < 1/(2Z). This 
corresponds to p > ph = 21.6 Z'^A g cm“^, Z and A being 
the charge and mass numbers of atomic nuclei. In this 
regime, Oe is considerably smaller than the iF-shell radius 
of an atom, and we have a fully ionized one-component 
plasma (OCP) of ions immersed in the “rigid” electron 
background. Then the pressure is 

P = Pe + Pi+Pc, (7) 


where Pg is the electron pressure. Pi is the pressure of ideal 
gas of ions, and Pq is the Coulomb term. 

The pressure of the ideal gas of electrons of any degen¬ 
eracy is given by (e.g., Landau & Lifshitz 1986), 


Pe = 


ksT 

Tr'^h^ 



1 -I- exp 


knT J 


p^Ap, 


( 8 ) 


where e = -^/m^cA + p'^c^ is the energy of an electron with 
a momentum p. The electron chemical potential p has to 
be determined from the equation 


1 p'^dp 

Jo 1+ exp[(e-p)/(fcBP)]' 


(9) 


The pressure contribution from the ideal gas of ions in 
Eq. d) is 


Pi = n-ik^T, 


( 10 ) 


where rii is the ion number density. 

The Coulomb correction accounts for electrostatic in¬ 
teractions of plasma particles. According to Hansen & 
Vieillefosse (1975), the results of Monte Carlo simulations 
of the OCP of ions immersed in a rigid electron back¬ 
ground can be fitted by 


Pc = Pi r3/2 


Ai 


-b 


Ao 


(Pi + r)i/2 B2 + r 


( 11 ) 


where Ai = -0.899962, Pi = 0.702482, A 2 = 0.274105, 
and P 2 = 1.319505. The fit error is smaller than 1% at 
r > 1, as confirmed by comparison with the recent Monte 
Carlo results of Stringfellow et al. (1990) and DeWitt et al. 
(1996). Equation ( 0 ) reproduces also the correct Debye- 
Hiickel limit at P <C 1. 

Equations (0)-0) cannot be used at rg ^ Z ^ and 
r ^ 1 because of the failure of the rigid electron back¬ 
ground approximation. They lead even to negative pres¬ 
sure at low temperatures and intermediate densities. Van 
Riper (1988) modified the Coulomb correction (0) in an 
ad hoc manner to avoid negative pressures. We use an¬ 
other approach described below in Sect. 2.2.4 


2.2.3. Low density regime 

In the low density regime, rg 3> 1 and T <C 1, the plasma 
is nearly ideal and may be partially ionized, depending 
on T and p. Theoretical description of partially ionized 
plasmas can be based either on the physical picture or 
on the chemical picture of the plasma (e.g., Ebeling et al. 
1977). In the chemical picture, the bound species (atoms, 
molecules, ions) are treated as elementary members of the 
thermodynamic ensemble, along with free electrons and 
nuclei. In the physical picture, nuclei and electrons (free 
and bound) are the only constituents of the ensemble. 
Both pictures can be thermodynamically self-consistent, 
though the chemical picture has limited microscopic con¬ 
sistency (see, e.g., Chabrier & Schatzman 1994, for re¬ 
view) . 

The most advanced results based on the physical pic¬ 
ture are those obtained in the OPAL project (Rogers et 
al. 1996). However, as an intrinsic limitation of the for¬ 
malism, they are restricted to low densities and/or high 
temperatures. We employ the OPAL EOS for partially 
ionized iron, which composes non-accreted envelopes. 

For the outermost layers of accreted envelopes, com¬ 
posed of H and He, we use the EOS derived by Saumon et 
al. (1995) (SCVH) within the framework of the chemical 
picture. The second-order thermodynamic quantities pro¬ 
vided by SCVH suffer from some thermodynamic incon¬ 
sistency in some regions of the phase diagram (see SCVH 
and Potekhin 1996, for a discussion), but the inconsistency 
is quite small at low densities. A comparison of the SCVH 
and OPAL EOSs for H and He reveals no discrepancies 
which could noticeably affect the temperature profiles in 
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a NS envelope under the conditions of interest. Further¬ 
more, the SCVH tables cover a larger p — T domain and 
provide also, along with the first- and second-order ther¬ 
modynamic quantities, the fractions of H and He atoms, 
H 2 molecules, and H+, He"*", He"*"’*' ions, which enables 
us to determine an average ion charge , requ ired for our 
calculations of the conductivities (Sect. 2.3.2 ). 


2.2.4. Intermediate density regime 

The intermediate density domain [Z~^ Ts 1) is most 
complicated because of the onset of recombination. Previ¬ 
ous studies of the thermal structure of the NS envelopes 
(GPE, Van Riper 1988) employed the model described by 
Eqs. with the free electron density Ue taken from 

the Los Alamos Astrophysical Opacity Library (Huebner 
et al. 1977, hereafter LAO). However, this approach fails 
at intermediate densities and relatively low temperatures 
(see, e.g.. Van Riper 1988), where it leads to unrealistic 
EOS owing to too large values of the Coulomb correction. 
In the afore-mentioned density region, the less dense and 
more strongly correlated electron gas is polarized by the 
external ionic field and eventually the electrons become 
bound by pressure recombination. Equation ( 0 ) becomes 
inadequate in this case, since it applies only to a fully 
ionized plasma, immersed in a rigid electron background. 

Fortunately, at temperatures of present interest 
(logT[K] ^ 4.7), the SCVH tables extend up to p > ph 
(~ 400 g cm“^ for He), i.e., to the completely ionized re¬ 
gion. Thus, they fully cover the intermediate density range 
for H and He. We use these tables up to the highest tab¬ 
ulated densities, and we use the high-density EOS given 
by Eqs. (0)-(0) for still higher p. 

The high-density domain for C and O starts from ph ~ 
10^ g cm“^. We do not need to consider C and O at lower 
densities in the present work. 

The case of iron is more difficult. The OPAL tables are 
bound from low-T high-p side approximately by the line 
P = POPAL = 10(T/10®K)^ g cm“^. In order to extend 
the EOS beyond this boundary, we adopt the approach of 
Fontaine et al. (1977), who circumvented a similar diffi¬ 
culty by interpolating pressure logarithm along isotherms 
over the intermediate density region. However, the conven¬ 
tional cubic spline interpolation of log P vs log p violates 
the condition {dP/dT)p > 0. Although this “normality 
condition” is not required by the thermodynamic consis¬ 
tency, it is supposed to hold at intermediate densities (see, 
e.g., Fontaine et al. 1977). 

In order to keep this condition, we introduce the fol¬ 
lowing modification of the interpolation procedure, (i) 
Pressure isotherms from logT = 4.5 to 8 are taken 
from the OPAL data at p < poPAL and calculated from 
Eqs. (|^)-([ll|) at p > Ph. (ii) The differences AlogP be¬ 
tween log P for neighbouring isotherms are interpolated 
across the gaps in logp between poPAL and ph by cubic 
splines, (hi) Starting with the isotherm logT(K) = 8.25, 



Fig. 1. Isotherms of effective charge for Fe obtained from the 
interpolated EOS at logr[K] = 4.5 (solid line), 4.75 (dashes), 
5 (solid), ... 8. 

for which Eqs. (|)-0) are sufficiently accurate, the pres¬ 
sure is consecutively calculated between poPAL and ph 
along the isotherms logr(K) = 8.0,7.75,...4.5, using 
the differences AlogP obtained at the preceding step: 
log P = log P' — A log P, where P' is the pressure at the 
preceding (hotter) isotherm. 

In order to verify this procedure, we have applied it to 
He and found that it gives much better agreement with 
the SCVH EOS than the straightforward logP-logp in¬ 
terpolation across the intermediate densities. 

2.2.5. Effective charge 

We will treat the electron heat conduction in the mean 
ion approximation, that is, we consider the plasma as 
composed of electrons and one ionic species with an ef¬ 
fective charge ZeS- In the case of H and He, we adopt 
the mean ion charge Z^s = (Z) provided by the SCVH 
tables, as discussed above. For carbon and oxygen, only 
the high-density, fully ionized regime is involved, for which 
ZeS = Z = 6 and 8, respectively. Finally, for iron we ad¬ 
just the effective number density of free electrons for 
the pressure Pe -I- Pi, calculated formally from Eqs. (|^- 
to reproduce the pressure P given by the accurate 
EOS at the same p and T. The mean ion charge, that 
corresponds to rie, is then assumed to be equal to Z^.^. 

Figure]^ shows the isotherms of Z^g corresponding 
to the interpolated Fe EOS. The decreasing parts of the 
curves reflect the electron recombination due to the reduc¬ 
tion of phase space per particle with increasing density. In 
this low-density region, our Z^fi coincides with {Z) given 
by the OPAL data (when available) within ~ 10%. The 
increase of Z^s at high densities corresponds to pressure 
ionization. The non-monotonic behaviour at intermediate 
densities reflects consecutive pressure destruction of dif¬ 
ferent electron shells of Fe atoms. 
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2.3. Opacities 

2.3.1. Radiative opacities 

In recent years, a significant progress has been made in cal¬ 
culations of the radiative opacities in dense atmospheres. 
We use the most recent OPAL opacity library (Rogers 
et al. 1996). For iron and hydrogen, it displays the ta¬ 
bles of the spectral opacities, which we convert into the 
Rosseland opacities RTrad using the standard frequency¬ 
averaging procedure (e.g., Schwarzschild 1958). For he¬ 
lium, we implement the Rosseland opacities readily avail¬ 
able in the OPAL library. 

For {p,T) values within the table range (p < popal), 
we have used bilinear interpolation of log iFrad between the 
table entries in logp and logT. For p > poPAL, we employ 
similar bilinear extrapolation based on the nearest pairs 
of table entries in logp and logT. The extrapolated values 
of IFrad have practically no effect on temperature profiles, 
because at high densities the heat transport is provided 
by electrons {Kc K^^d). 

2.3.2. Conductive opacities 

Previous studies of the temperature profiles in NS en¬ 
velopes have mainly used the electron thermal conduc¬ 
tivity of completely ionized, strongly degenerate plasma. 
In particular, one has often taken the conductive opaci¬ 
ties from Urpin & Yakovlev (1980) and Yakovlev & Urpin 
(1980) (hereafter YU). However, in a cold enough NS, the 
thermal conductivity of non-degenerate or partly degener¬ 
ate electrons may become important. A comparison with 
the tabular data of Hubbard & Lampe (1969) reveals that 
a straightforward extrapolation of the YU formulae from 
their validity domain (fully ionized, degenerate plasma) to 
the case of non-degenerate matter may lead to an underes¬ 
timation of the electron thermal conductivity by orders of 
magnitude. As will be shown in Sect. this may strongly 
affect the temperature profiles in the NS envelopes at 
Teff ^ 105 K. 

In order to determine the conductive opacities for any 
degeneracy, we employ the numerical code developed re¬ 
cently by Potekhin & Yakovlev (1996) (hereafter PY) for 
calculating the electron transport coefficients along mag¬ 
netic fields in the magnetized NS envelopes. The code per¬ 
forms numerical thermal averaging of the effective energy- 
dependent electron relaxation time with the Fermi-Dirac 
distribution at any degeneracy. This enables us to consider 
not only strongly degenerate, but also mildly degenerate 
matter. In addition, PY have improved the approach of 
YU by a more accurate calculation of the Coulomb loga¬ 
rithm for the electron-ion scattering in the liquid phase of 
NS matter and by incorporating the Debye-Waller reduc¬ 
tion factor (e.g., Itoh et al. 1984) for the electron-phonon 
scattering in the solid phase. 

In the present article, we apply the above code for the 
particular case of zero magnetic field. For this purpose. 


further modifications have been made. First, in addition 
to the thermal conduction produced by the electron-ion 
scattering considered by PY, we have incorporated also 
the contribution from the electron-electron scattering in 
the liquid regime (see Sect. [A.lD . This process can be im¬ 
portant for light elements (H, He) in a weakly degenerate 
matter. Second, we have improved the Coulomb logarithm 
which determines the collision rate o f strongly degenerate 
electrons with ions (see Sect. 2.3.3| ). Third, in the solid 
regime, we have used the new expression for the elec¬ 
tron thermal conductivity, derived recently by Baiko & 
Yakovlev (1995). Note that at temperatures much lower 
than the melting temperatures, the Coulomb scattering 
of electrons by charged impurities may become important 
(e.g., YU). It depends on the impurity charge and, most 
important, on its number density Uimp which is generally 
unknown in the NS envelopes. We will present the results 
obtained for pure crystals (nimp ^ 0). Our calculations 
including ^^Ni impurities in iron envelopes show that the 
electron-impurity scattering may affect the temperature 
profiles in the crusts of cooling NSs, whenever the im¬ 
purity concentration is sufficiently high, nimp ^ 0.01 ni. 
Finally, in the non-degenerate layers, we have performed 
the thermal averaging introducing modifications into the 
effective width of the Landau levefs which was used by 
PY to account for coflisional and inefastic broadenings of 
the Landau levefs in strongly degenerate matter. The ef¬ 
fective width proposed by PY becomes inadequate in the 
non-degenerate regime, and we have switched it off by 
multiplying by the factor (1 -|- 10T/Tf)“^. 

A comparison of our thermal conductivities with the 
tables of Hubbard & Lampe (1969) (at p ^ lO® g cm”^^ 
where the tables are adequate) shows maximum discrep¬ 
ancy up to 30% for carbon and up to 70% for lighter 
elements in the non-degenerate regime, and still smaller 
discrepancies in mildly and strongly degenerate matter. 
More a ccurat e values of the Coulomb logarithm obtained 
in Sect. 2.3.3 modify Kc by about 5-30%, which has only 
a little effect on the thermal structure of the envelope (see 
Sect. H). 


2.3.3. Coulomb logarithm 

The thermal conductivity k of degenerate electrons {T < 
Tp) in a liquid or a gas of ions with mass and charge 
numbers A and Z can be written as 


TT^/CgTrie 


( 12 ) 


where m* = WeVl + Vc = + Vee is the effec¬ 

tive electron collision frequency, t'ei is determined by the 
Coulomb electron-ion collisions, and Vge is determined by 
the electron-electron collisions. New calculations and fits 
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of Vee are presented in Sect. A.l. In this section, we exam¬ 
ine t'ei- It can be expressed as (e.g., YU) 


^ei 




(13) 


where L is the Coulomb logarithm to be evaluated. It is a 
slowly varying function of p and T. 

First consider a strongly coupled Coulomb liquid of 
ions, at 1 ^ r Fm, where Fm is the critical value of 
the ion coupling parameter F at which a Coulomb crys¬ 
tal melts (Fni = 172, for a classical ion system, Nagara et 
al. 1987; the deviations from the classical melting curve 
due to the quantum corrections turn out to be small at 
p < 10^° g cm“^ for carbon and heavier ions - see Chabrier 
1993). In this case, one usually assumes that the screen¬ 
ing of the electron-ion Coulomb interaction is static and 
produced by the ion-ion correlations. The screening is de¬ 
scribed by the ion-ion structure factor S (q). Extensive cal¬ 
culations of L in the Born approximation for different ion 
species, Z and A, were performed by Itoh et al. (1983). 
The authors used the structure factors obtained in the 
approximation of rigid electron background. They fitted 
their results by complicated expressions. Later Yakovlev 
(1987) evaluated non-Born corrections to L, using the 
weak screening approximation. 

Now we recalculate the Coulomb logarithm for 
strongly degenerate electrons at 1 ^ F ^ Fm from the 
expression 


L = 


q^Siq)F^iq)R{q) 

3 


1 - 




(14) 


where (3 = xj\J\ -I- kp = pp/Ti, Tiq is a momentum 
transfer in an ei collision event, e{q) is the static longitu¬ 
dinal dielectric function of electrons, F{q) is the nuclear 
form factor to allow for finite nuclear size, and R{q) is 
the non-Born correction factor specified below. The di¬ 
electric function has been calculated for the electron gas at 
any degeneracy. However, we have verified that, for all the 
cases in the present study, this dielectric function yields 
the same results as the function obtained in the limit of 
strong electron degeneracy for the relativistic electron gas 
(Jancovici 1962). Thus the bulk of the calculations has 
been done with this latter function. We use the nuclear 
form factor F{q) appropriate to a uniform proton core of 
an atomic nucleus, with the same core radii as in Itoh et 
al. (1983). The finite sizes of the nuclei appear to have 
almost no effect on L, for the conditions of study. Never¬ 
theless we have kept F{q) in our calculations. Following 
Yakovlev (1987), the non-Born factor has been taken as 
R{q) = (T(q)/(TB(g), where a(q) is the exact differential 
Coulomb scattering cross section for a momentum trans¬ 
fer g, and <Jp,{q) is the cross section in the Born approxi¬ 
mation. The exact cross sections are taken from Doggett 
& Spencer (1956). 


In our calculations from Eq. (0 , we have used the 
ion structure factors evaluated for a responsive electron 
background including the local field corrections between 
electrons (Chabrier 1990). Our calculations extend those 
of Itoh et al. (1983) by including the non-Born correc¬ 
tions and the effects of electron response on the ion struc¬ 
ture factors. On the other hand, we improve the results 
of Yakovlev (1987) since we evaluate the non-Born correc¬ 
tions beyond the weak screening approximation. 

We have evaluated L from Eq. & for chemical ele¬ 
ments with 1 < Z < 26 and for a dense grid of 1 < F < Fm 
and 10^ g cm“^ F p < 10^° g cm“^ provided T < Tp. 
In the case of hydrogen, we have restricted ourselves to 
p < 10® g cm“® since hydrogen should be completely burnt 
at higher densities. For light elements, the highest tem¬ 
peratures included in these data correspond to F = 1 and 
appear to be much below Tp. Accordingly, there exists a 
large temperature interval in the electron degeneracy do¬ 
main where F < 1, the ion coupling ceases to be strong, the 
ion screening cannot be treated in terms of ion-ion corre¬ 
lations, and Eq. (Q) is invalid. In principle, the Coulomb 
logarithm for mildly-coupled ion plasma (F ^ 1) can be 
evaluated using the formalism developed by Boerker et al. 
(1982). In order to simplify our analysis, we will restrict 
ourselves to the weak-coupling case (F 1), in which the 
ions constitute a Boltzmann gas, the weak ion screening 
is dynamical and can be described by the dielectric func¬ 
tion formalism, while the electron screening remains static 
(e.g., Lampe 1968). At F <C 1 and T < Tp, we propose 
the expression 


- + \ nZaP, 




(15) 


where b, = l/{2kprp,i) = l/y/^ is the ion screening pa¬ 
rameter in the weak coupling regime, rp,-, is the Debye 
radius of ions, p = be/h,, be = kTp/(2kp) = \Jaj (tt/?) is 
the electron screening parameter, fcxF being the inverse 
screening length of a charge by degenerate plasma elec¬ 
trons and a the fine-structure constant. The first three 
terms represent the Coulomb logarithm obtained for de¬ 
generate electrons and weakly coupled ions with the non- 
relativistic ei scattering cross section. The expression pre¬ 
sented is valid for ?; < 1; it can be derived using the for¬ 
malism of Williams & DeWitt (1969) to account for the 
dynamic ion screening and to integrate the ei collision rate 
over velocities of ions. The fourth term is the relativistic 
correction in the weak-coupling limit (e.g., YU). The fifth 
term is the second-Born correction in the weak-coupling 
limit (Yakovlev 1987). Thus, we have supplemented our 
table of the Coulomb logarithms calculated for F > 1 
(see above) with new values evaluated from Eq. (15) for 
F < 0.25 and T < Tp. Analytic fits to the full set of about 
31,000 values of L are presented in Sect. |A.2| . The tables 
of L are freely available in the electronic form. 
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Fig. 2. Coulomb logarithm L vs density for fully ionized plas¬ 
mas at r = 10. Left panel: Fe and C plasmas; dash-an d-dot s 
— present calculations, solid line — analytic fits (Sect. [Y^ ), 
dashes — Born approximation (Itoh et al. 1983). Right panel: 
H and He plasmas; dots — present calculations, solid line — 
analytic fits, dash-and-dots — ion-structure factors in the rigid 
electron background (Itoh et al. 1983). 


For illustration, in Fig. || (left panel) we compare the 
Coulomb logarithms in fully ionized Fe and C plasmas 
with those calculated in the Born approximation (Itoh et 
al. 1983). The non-Born corrections enhance the small- 
angle Coulomb scattering cross section, and increase L. 
The increase is especially pronounced for the heavier el¬ 
ement, Fe (exceeds 20%), at densities p ^ 10® g cm“®, 
where the electrons are relativistic. The non-Born correc¬ 
tions become less important with decreasing Z, and they 
are negligible for light elements (H, He). 

The inclusion of the electron response in the calcula¬ 
tions of the ionic structure factors affects the Coulomb 
logarithm L for H and He at F ^ 1, as demonstrated on 
the right panel of Fig. ^ One can observe an increase of 
L up to 15% for H at p = 10^ g cm“^ and F = 10 with 
respect to the case of rigid electron background. Actu¬ 
ally, the increase grows with F and reaches about 40% at 
F = 160. The effect is weaker for He, and small for heav¬ 
ier elements (in the adopted parameter range). Thus, for 
the elements heavier than He, we have used the structure 
factors calculated in the rigid background approximation. 

Note that calculated values of L determine also the 
electric conductivity of degenerate matter in a liquid or 
gaseous ionic plasma, a = e^ne/(TO*r'ei). 


3. Results and discussion 

3.1. Thermal structure of non-accreted envelopes 

In this section, we consider a non-accreted NS envelope 
composed solely of iron, which may be at various ioniza- 



log r, (K) 


Fig. 3. Boundary conditions for the thermal structure equa¬ 
tion (|^. Surface density is obtained from Eq. (^ as a func¬ 
tion of surface temperature, for gravities pi4=0.4, 1, 2.43, 6. 
Dash-dotted lines correspond to a non-accreted matter (Fe 
opacities), and solid lines to an accreted matter (H opacities). 


S) 
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tion stages. We use the OPAL iron EOS where available, 
the ideal-gas EOS with the Coulomb correction ( 0 ) at 
high densities, and the interpolated EOS at intermediate 
densities as described in Sect, p.2.4 . The conduction prop¬ 
erties are described in the mean ion approximation with 
the effective ion charge determined consistently with 
the EOS (Sect. 12.2.51) . 

Figure ^ presents density dependence of the temper¬ 
ature in the NS envelope at various Tefj. The integration 
of Eq. (|^) has been started with the surface density p, 
shown in Fig. and terminated at pb = 10^® g cm 
The value gu = 2.43 chosen in Fig. ^ corresponds to a 
NS with a radius i? = 10 km and a mass M = 1.4 M©. 
Circles on the curves are the points of equal radiative and 
conductive opacities. The thermal flux is mainly carried 
by radiation at lower densities and by electrons at higher 
densities. We also show the electron degeneracy curve and 
the melting curve. 

In a wide range of Tea, some layers (Fig. turn out 
to be convectively unstable (Zavlin et al. 1996, Rajagopal 
& Romani 1996). In these layers, the energy transport is 
dominated by convection rather than by electron or radia¬ 
tive conduction. We describe the convective energy flux in 
the adiabatic approximation (Schwarzschild 1958), which 
assumes that the convective heat transport is much more 
efficient than other mechanisms. Accordingly, in the con¬ 
vective zones, we replace the temperature gradient given 
by Eq. (|^ by the adiabatic gradient provided by the EOS 
tables. 

In order to check the effect of convection, we have per¬ 
formed also calculations with “frozen” convection (i.e., ne¬ 
glecting the convection effect). This is the extreme case 
opposite to the adiabatic one. In this approximation, we 
obtain (Fig. ||) slightly higher temperatures inside the con¬ 
vective part of the atmosphere. The atmospheric temper- 
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Fig. 4. Temperature profiles in a non-accreted NS (heavy solid 
lines) compared with various approximations: after Yakovlev & 
Urpin (1980) (YU) by extrapolating the conductive opacities 
to T > Tf; after Potekhin & Yakovlev (1996) (PY); and with 
neglected convection (thin solid lines). The curves are labeled 
by the values of Tefr/10® K. Circles show the points, where the 
radiative opacity equals the conductive one; heavy dots show 
the melting curve P = 172; long dashes display the degeneracy 
curve, T = Tp. 


ature profiles derived by Zavlin et al. (1996), who solved 
numerically the radiative transfer equation at moderate 
optical depths and described the convection using the 
mixing-length theory, lie between our two extremes. In 
deeper layers, the two extreme profiles tend to merge, as 
clearly seen in Fig. |^. This merging occurs because the fac¬ 
tor K/T'^ on the r.h.s. of Eq. ( 0 ) rapidly decreases, with 
increasing T, and thus ensures the isothermal profile at 
higher T behind the convective zone (at fixed Teff). Thus, 
the thermal structure of the NS envelope at p ^ 10 g cm“^ 
is practically unaffected by the convection. 


The dot-dashed curves in Fig. |4| are obtained with Kc 
calculated according to PY (Sect. ]23^ ). These curves al¬ 
most coincide with the solid lines, calculated more accu¬ 
rately. This confirms the validity of the PY approxima¬ 
tions for the non-accreted envelopes. 


If Teff ^ 2 X 10® K, the electron conduction becomes 
important not only in degenerate matter, but also in a 
non-degenerate region (Fig. Then the straight-line log¬ 


Fig. 5. Temperature increase through non-accreted NS en¬ 
velopes with different values of the surface gravity as compared 

to Eq. ( 0 ). 


arithmic extrapolation of the YU conductive opacities into 
the non-degenerate regions leads to a significant overesti¬ 
mation of the internal temperature (dotted curves). 

This effect is seen also in Fig. |^, that shows 
log(Tb/reff) as a function of log Teg for four values of the 
surface gravity. If logTs [K] > 5.5, the well-known simple 
fit of GPE, 

rb6 = 128.8 (r4/5i4)° «®, (16) 

is rather accurate, but at lower Tg the deviations from 
this fit are quite pronounced. Nevertheless, the scaling law 
7b = f{Ts/g) holds with a good accuracy in the whole 
temperature-gravity range presented in this figure, except 
for the lowest Tb values. The fit which reproduces these 
new features is presented in Sect. |A.3| . 

3.2. Thermal structure of accreted envelopes 

The chemical composition of an accreted envelope de¬ 
pends on the composition of the infalling matter and 
on the nuclear transformations (thermo- or pycnonuclear 
burning, beta-captures) which occur while the matter 
sinks within the NS. The structure and stability of the ac¬ 
creted envelopes have been considered in number of works 
devoted to bursting NSs (see, e.g., Ayasli & Joss 1982, 
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log density (g/cc) 


Fig. 6. Temperature profiles in an accreted NS (solid lines) 
as compared to those in NS envelopes composed of pure iron 
(dot-dashed lines) and of pure helium (dashed line). The curves 
are labeled by Teff/IO^K; full circles separate the regions of 
radiation and electron conduction; asterisks indicate the H/He 
(left) and He/C (right) interfaces. 


Paczynski 1983, Ergma 1986, Miralda-Escude et al. 1990, 
Blaes et al. 1992, and references therein). However, most 
of the work focused either on high accretion rates or on 
slowly accreting old NSs whose internal temperature is 
solely determined by the accretion. On the contrary, we 
are mainly interested in not too old, cooling NSs with 
thin accreted envelopes, in which the heat release due to 
the nuclear transformations has no effect on the thermal 
structure. 

Let us assume that an accreted envelope consists of 
shells with different chemical compositions. Eor an accre¬ 
tion of H/He, it is reasonable to expect that the outer¬ 
most accreted layers are built of pure H, owing to the 
strong gravitational stratification (Alcock & Illarionov 
1980), whereas He may exist in the deeper layers. The ac¬ 
cretion is accompanied by the nuclear transformations of 
H/He into heavier elements. The details of these transfor¬ 
mations are still uncertain, and we consider several models 
of accreted envelopes. 

As a first example, we take the typical temperatures 
and densities of the nuclear burning from Ergma (1986) 


and assume the burning to be non-explosive. If Tgff is high 
and the accreted hydrogen layer reaches the region with 
T ~ 4 X 10^ K, then H burns efficiently into He in the 
thermonuclear regime. If temperature is lower and the 
hydrogen layer reaches the density p ~ 10^ g cm“^, a 
pycnonuclear burning of H into He proceeds. At higher 
p ~ 10® g cm“^ and/or T ^ 10® K, helium burns into 
carbon (Schramm et al. 1992). The pycnonuclear carbon 
burning occurs at p ^ 10^® g cm“® (e.g., Yakovlev 1994), 
which is already inside the isothermal region of the enve¬ 
lope. 

Figure ^ displays the thermal structure of a fully ac¬ 
creted envelope (the accreted mass AM ~ 10“^M fills 
the layers up to p ~ pb) of a NS with the surface gravity 
Ps = 2.43 X 10^"* cm s“®. The outer, intermediate, and 
inner shells of this envelope, separated by asterisks, are 
composed of H, He, and C, respectively. 

One can observe significant differences from the ther¬ 
mal structure of a non-accreted Fe envelope, as explained 
below. For a not too cold NS (T/ff ^ 10® K), the main 
temperature gradient occurs in a layer of degenerate elec¬ 
tron gas with ions in the liquid state, where the ther¬ 
mal conduction is mostly provided by the electrons due to 
the electron-ion collisions. The heavier the element, the 
smaller the thermal conductivity, and the steeper is the 
temperature growth inside the star. However, with de¬ 
creasing Teff, the width of the heat-insulating degenerate 
layer becomes smaller, and the effect is less pronounced. 
In a cooler NS {Tes ^ 10® K), the main temperature gra¬ 
dient shifts into the NS atmosphere (to the optical depths 
^ 1). For heavier elements, the atmospheric layers appear 
generally denser (at the same T^g, see Fig. ||). Then the 
internal temperature gradient is weaker and the temper¬ 
ature grows slower inside the star. Thus, a not too cold 
accreted envelope is more heat-transparent than the iron 
envelope, whereas a cold accreted envelope is less heat- 
transparent. As seen from Fig. |^, the effective surface tem¬ 
perature which separates these two regimes is almost in¬ 
dependent of the surface gravity. 

Figure 1^ demonstrates the effect of the improvements 
in the conductive opacities. The present results are com¬ 
pared with those in which the electron thermal conduc¬ 
tivity is calculated with the PY code. Although the differ¬ 
ence is not too large, it is higher than for the non-accreted 
Fe envelopes (Fig. The difference comes mostly from 
the electron-electron collisions (neglected in PY). Dots in 
Fig. 0 show the results obtained with the PY code sup¬ 
plemented by the contribution from the electron-electron 
collisions to the thermal conduction (Sect. A.l). The dif¬ 
ference from the accurate results practically disappears. 

Now let us analyze the sensitivity of the temperature 
profiles to the positions of the H/He and He/C interfaces. 
These interfaces, indicated by asterisks in Fig. |[ are not 
well theoretically established. For example, the densities of 
nuclear H and He burning presented by Iben (1974) (see 
also Kippenhahn & Weigert 1990) are much lower than 
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log p [g cm-:=] log T^„ [K] 


Fig. 7. Temperature increase through accreted NS envelopes 
for selected values of the surface gravity calculated using the 
full code (solid lines), PY-code (dot-dash lines), and PY-code 
supplemented by thermal conduction due to electron-electron 
collisions (dots). 


those given by Ergma (1986). We have significantly shifted 
the interfaces according to the results of Iben (1974), but 
the Teff-Tb relationship remains practically unchanged. 
We have also calculated the temperature profiles for a 
purely helium NS envelope. These profiles are slightly dif¬ 
ferent from those in a burning accreted matter described 
above. However, the Teff-Tb relationships remain practi¬ 
cally the same. 

Recently the thermal structure of old, slowly accreting 
NSs has been examined in several articles (see Miralda- 
Escude et al. 1990, Blaes et al. 1992, and references 
therein). In these studies, the internal stellar tempera¬ 
ture is assumed to be determined by the accretion rate. 
The composition of the accreted envelopes appears to be 
similar to the ones considered above except for the in¬ 
nermost layers. For instance, Miralda-Escude et al. (1990) 
assume that He transforms directly into Fe while Blaes et 
al. (1992) state that the burning of He into C at p ^ 10® 
g cm“^ is followed by the carbon transformation through 
the rapid ^®C(a,7)^®0 process. Accordingly, we have re¬ 
calculated the temperature profiles replacing by ®®Fe 
or and found only very small differences. The insensi¬ 
tivity of the temperature profiles to the composition of the 
innermost layers of the envelope is quite natural, because 
these layers are nearly isothermal. 

Summarizing the above discussion we can conclude 
that the Teff-Tb relationship is remarkably insensitive to 
the parameters of nuclear burning and to the chemical 
composition of accreted matter (H and/or He) but differs 
from the one in a non-accreted NS. This conclusion will 
not change if H burning is explosive (see, e.g., Paczynski 
1983, Miralda-Escude et al. 1990). During an outburst. 



Fig. 8. Temperature increase through partly accreted NS 
crusts. 


the H-shell burns into He, without any noticeable effect 
on the Teff-Tb relation. According to Miralda-Escude et 
al. (1990) helium burning in slowly accreting NSs is un¬ 
stable; it produces nuclear outbursts with burning of all 
accreted material. Nevertheless, the helium layer should 
reach a very high density, p ^ 10® g cm“^ to trigger the 
instability in a not too hot NS. If the same is true for the 
conditions of our interest, our results are valid during a 
long period of time until the instability is reached. How¬ 
ever, according to Paczynski (1983), the presence of the 
heat outflow from a not too cold, cooling NS essentially 
stabilizes nuclear burning in accreted envelope, and the 
He burning is mostly non-explosive for the parameters we 
are interested in. 


Figure^ shows log(Tb/Teff) as a function of logTeff 
for various amounts of accreted matter. The dot-dashed 
line represents log(Tb/Teff) for a non-accreted envelope 
from Fig. ||. This function is strongly affected by the ac¬ 
cretion. Even a thin hydrogen (or He) mantle of mass 
AM = 10“^®M, which extends only to p ^ 10^ g cm“^, 
modifies strongly the Tb-Teff relation. Analytical fits of 
Teff as a function of the internal temperature Tb, surface 
gravity g, and accreted mass AM are given in Sect. A.3. 
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3.3. Effect of accreted envelope on cooling 

Consider briefly the effects of the possible presence of ac¬ 
creted matter in the NS envelopes on the thermal history 
of NSs. For this purpose, we simulate the cooling of NSs 
using two NS models. We use the cooling code of Gnedin 
& Yakovlev (1993) based on the approximation of isother¬ 
mal NS interior. The code follows accurately the cooling 
of a not too young NS whose thermal relaxation is already 
over {t ^ 10^ yrs, see Nomoto & Tsuruta 1987, Umeda et 
al. 1994). We adopt a moderately stiff EOS of Prakash et 
al. (1988) to describe matter in the NS cores. This matter 
is assumed to consist of neutrons, protons, and electrons 
(no exotic cooling agents such as quarks or kaon conden¬ 
sates). The maximum NS mass, for this EOS, is about 
1.7 Mq. We examine two NS models: (1) M = 1.30 Mq, 
R = 11.72 km, the central density Pc = 1-12 x 10^® g 
cm-3; (2) M = 1.44 Mq, R = 11.35 km, = 1.37 x lO^^ 
g cm“^. 

The first model is a typical example of the standard 
cooling. The central density is insufficient to switch on 
the most powerful neutrino cooling due to the direct Urea 
process (Lattimer et al. 1991). In this case the main neu¬ 
trino energy loss mechanisms in the NS cores are the mod¬ 
ified Urea reactions (neutron and proton branches) and 
the nucleon-nucleon bremsstrahlung (Friman & Maxwell 
1979, Yakovlev & Levenflsh 1995). 

Our second NS model is an example of the rapid neu¬ 
trino cooling. In this model, the central density slightly 
exceeds the threshold density pa = 1.30 x 10^® g cm“^ 
of the direct Urea process (for the given EOS). Then the 
NS has a small central kernel (of radius 2.32 km and mass 
0.035 Mq) where the direct Urea process operates, and the 
neutrino luminosity exceeds the standard one by several 
orders of magnitude. 

All the neutrino energy loss rates included in the cal¬ 
culations are described by Levenflsh & Yakovlev (1996). 
To simplify our analysis, we assume that the NS cores are 
non-superfluid, and there is no internal reheating (e.g., en¬ 
ergy dissipation due to differential rotation). The results 
are presented in Fig. in the form of the cooling curves 
- effective temperatures as measured by a distant ob¬ 
server (see Sect. |2.lD versus NS age t. 

Figure ^presents the cooling curves under the assump¬ 
tion that all the accreted material is accumulated at the 
NS surface during the NS birth. The accretion can oc¬ 
cur, for instance, at the post-supernova stage (e.g., Cheva¬ 
lier 1989, 1996, Brown & Weingartner 1994) especially if 
the pulsar activity is suppressed initially by burying the 
pulsar magnetic held under the fallen-back matter (Mus- 
limov & Page 1995). We show the standard and rapid neu¬ 
trino cooling for different amount of accreted mass AM. 
The fraction of accreted mass AM/M varies from 0 (non- 
accreted NSs) to ~ 10“^ (fully accreted NS envelope). It 
is evident that further increase of AM does not affect the 
cooling. 



Fig. 9. Standard and rapid cooling of NSs with non-accreted 
(Fe) and partly and fully accreted envelopes, compared with 
observations of thermal radiation from 6 pulsars (see text). 
The accreted mass is assumed to be accumulated during the 
first 100 yrs; numbers along the curves denote the values of 
log AM/M. 


The change of slopes of the cooling curves at t = 10^- 
10® yrs reflects the change of the cooling regime. Initially, a 
NS is at the neutrino cooling stage. It cools mainly via the 
neutrino emission; the internal stellar temperature is ruled 
by this emission and is thus independent of the thermal in¬ 
sulation of the NS envelope. The surface photon emission 
is determined by the relation. Since the accreted 

envelopes of not too cold NSs are more thermally trans¬ 
parent (Sect. ^), the surface temperature of an accreted 
NS is noticeably higher than that of a non-accreted NS. 
One can see that even a very small fraction of accreted 
matter, such as AM/M ~ 10“^®, may change appreciably 
the thermal history of the star. The colder the star, the 
smaller is the fraction of accreted material which yields the 
same cooling curve as the fully accreted NS envelope. The 
effect is naturally more pronounced for the rapid cooling 
(for colder NSs). 

If t ^ 10®-10® yrs, the neutrino luminosity decreases, 
and the NS cools mainly via the surface emission of pho¬ 
tons (the photon cooling stage). Since the envelopes of the 
accreting (and not too cold) NSs have lower thermal insu¬ 
lation, these stars cool now faster than the non-accreted 
ones (Fig. ^). 
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The change of the thermal history of a NS by the ac¬ 
creted material should be taken into account in theoret¬ 
ical interpretation of observational data on thermal ra¬ 
diation from NSs. For illustration, in Fig. we present 
the effective temperatures with Icr confidence intervals of 
the pulsars 0833-45 Vela (Ogelman et al. 1993, Ogelman 
1995), 0630-1-18 Geminga (Halpern & Ruderman 1993), 
0656-1-14 and 1055-52 (Greiveldinger et al. 1996), RXJ 
0002-1-6246 (Hailey & Craig 1995), and the 3a upper limit 
for the Crab pulsar (PSR 0531-1-21, Becker & Aschenbach 
1995). These data have been deduced from the fits of spec¬ 
tral observations of the ROSAT X-ray observatory. For 
Geminga, PSR 1055-52, and PSR 0656-1-14, we plot the 
values of Tes obtained for “soft” components of the two- 
and three-component X-ray spectra models, since it is the 
soft component that is associated with the thermal flux 
from the stellar interior. For Vela, we extend downward 
the confidence interval obtained by Ogelman et al. (1993) 
according to the more recent spectral fitting of Ogelman 
(1995), that gives a more plausible NS radius than the 
previous fit. 

The age of the Crab pulsar in Fig. ^ is its true age 
known from the chronicles, the age of Vela is its possible 
true age determined recently by Lyne et al. (1996), the age 
of RXJ 0002-1-6246 is an estimate for the related supernova 
remnant G 117.7-1-06, while for the three other pulsars we 
adopt their characteristic ages. 

We plot the values of Tea derived by the authors by 
fitting the assumed blackbody thermal spectrum to the 
measured X-ray spectrum. The alternative often used is to 
derive Tea from the apparent luminosity in the whole spec¬ 
tral range observed. However, this latter way is less cer¬ 
tain, because the inferred temperature strongly depends 
on the assumed NS distance (usually poorly known) and 
model-dependent radius. 

Note that the blackbody fits do not take into account 
the influence of the magnetized atmosphere. The fits for 
PSR 0656-1-14 (Anderson et al. 1993), Geminga (Meyer 
et al. 1994), and Vela (Page et al. 1996), which employ 
the models of magnetized atmospheres, have resulted in 
a few times lower Teff than the blackbody models. Our 
justification of using here the blackbody data is that the 
so far developed models of magnetized atmospheres do 
not take proper account of absorption by neutral atoms. 
The neutral fraction can be significant at T ^ 10® —10® K, 
because the magnetic fields ~ 10^^ G increase the ground- 
state binding energy by about one order of magnitude. 
The opacities of hydrogen in such strong fields are still 
under investigation (see Pavlov & Potekhin (1995) and 
references therein). 

The magnetic fields affect also the thermal insula¬ 
tion of the NS envelopes and, therefore, the cooling (e.g., 
Yakovlev & Kaminker 1994). We do not take this effect 
into account in the present article. As suggested by Hern- 
quist (1985) and shown recently by Shibanov & Yakovlev 
(1996), the dipole surface magnetic fields ^ 3 x 10^^ G af¬ 


fect the cooling much more weakly than predicted by the 
studies of simplified magnetic field geometries (e.g.. Van 
Riper 1988, 1991). 

As seen from Fig. for the pulsars Vela, Geminga, and 
RXJ 0002-1-6246, the presence of accreted matter would 
simplify the explanation of the cooling by the standard 
neutrino emission without invoking superfluidity of NS 
matter, as favoured in previous work (Page 1994, Lev- 
enfish & Yakovlev 1996). Interestingly, the accretion rate 
of Vela estimated by Morley (1996), multiplied by its age, 
yields log AM/Mq « —9.5, that just fits the middle of the 
error bar in Fig. For PSR 0656-1-14, the scenarios with 
and without the accretion are equally consistent with the 
result of the htting shown if Fig. |^, especially if we assume 
that its braking index is somewhat lower than 3 (as hap¬ 
pens for all 4 pulsars with known braking indices, see Lyne 
et al. 1996) and therefore its true age is somewhat higher 
than its characteristic age. The relatively high tempera¬ 
ture of PSR 1055-52 cannot be explained by the consid¬ 
ered models of the cooling, that probably suggests some 
reheating process for this relatively old pulsar. 

We have used, for illustration, the standard NS mod¬ 
els and simplified blackbody fits to the observed spectra. 
More sophisticated models and fits can also gain from 
the assumption of accreted envelopes. For instance, Page 
(1996) considers the cooling of NSs with kaon condensed 
cores and reheating in the outer crusts, using our models 
of the envelopes. 


4. Conclusions 


We have considered thermal structure and evolution of 
NSs whose envelopes are composed of non-accreted or ac¬ 
creted matter. We have used new, state-of-the-art calcu¬ 
lations of EOS and opacities of NS envelopes, descri bed in 
Sect . 2. In pa rticular, we have recalculated (Sects. |2.3.3 


A.l , and |A.2|) the electron-electron and electron-ion colli¬ 


sion frequencies, that determine the electron thermal con¬ 
ductivity, for a wide range of densities and temperatures 
of a degenerate electron gas and ionic liquid plasmas. 

Usi ng t his new physics input, we have calculated 
(Sects. the temperature profiles in the envelopes 

of non-accreted and partly accreted NSs and obtained the 
relationships between the internal and effective temper¬ 
atures of NSs, Tb and Teff. These relationships are im¬ 
portant for simulating the NS cooli ng; t hey are fitted by 
simple analytic expressions in Sect. AM. They appear to 
be very sensitive to the presence of accreted matter in the 
NS envelope; even a small amount of the accreted matter, 
AM ~ 10“^® Mq, can reduce substantially the thermal 
insulation of the envelopes. For a non-accreting NS, our re¬ 
lationship between Tb and Teff extends the well-known re¬ 
sult of GPE to lower temperatures (down to Teff ~ 50 000 
K). 


In Sect. 3.3, we have examined briefly the effect of the 
possible presence of accreted matter on the NS cooling. 
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We show that the accreted matter may increase the sur¬ 
face temperature (photon thermal luminosity) at the neu¬ 
trino cooling stage, and decrease them at the subsequent 
photon cooling stage, as compared to the NSs without ac¬ 
creted envelopes. We have shown that these results can be 
important for a proper interpretation of observed thermal 
radiation from NSs. In particular, the presence of accreted 
matter facilitates the explanation of recent observational 
results concerning the pulsars Vela and Geminga, and RXJ 
0002-1-6246, in the framework of the standard neutrino 
emission model (without exotic matter, superfluidity, or 
direct Urea processes). 

In this paper, we have neglected effects of magnetic 
fields on the EOS and the thermal conduction of mat¬ 
ter, which can be significant (e.g., Yakovlev & Kaminker 
1994). They do deserve further studies using improved 
EOS and thermal conductivities of magnetized NS enve¬ 
lope (e.g., PY) and improved radiative opacities of mag¬ 
netized NS atmosphere (e.g., Pavlov & Potekhin 1995), 
with allowance for the possible presence of light elements 
in the surface layers. 

Pinally, it is worth noting that the physics input used 
in the present calculations can be applied to a variety 
of other astrophysical problems concerning dense stellar 
matter, e.g. the thermal structure and bursting activity 
of X-ray bursters (see, e.g., Miralda-Escude et al. 1990 
and references therein) and the cooling of white dwarfs 
(Segretain et al. 1994). 

Acknowledgements. We are grateful to H. DeWitt for pointing 
out the article of Boerker et al. (1982), to Dany Page and Yu.A. 
Shibanov for fruitful discussions, and to D. Page for reading 
the manuscript and making the most useful critical remarks. 
This work was supported in part by the RBRF grant No. 96- 
02-16870a, INTAS grant No. 94-3834, and DFG-RBRF grant 
No. 96-02-00177G. AYP gratefully acknowledges the hospital¬ 
ity of the theoretical astrophysics group of the Ecole Normale 
Superieure de Lyon. 


Appendix: fitting formulae 


(1976) to higher temperatures, T < Tp. In the approxi¬ 
mation of static electron screening of the Coulomb inter¬ 
action, Urpin & Yakovlev (1980) obtained 


^ee 


{kBTf 
27r3 hmlc^ 


J{y), 


(A2) 


where y = V^Tpe/T, a = e^/{hc), /3 = x/Vl + and 
be = a/^irP) (see Eq. (15)). Now it is sufficient to calcu¬ 
late the function J{y), presented by Urpin and Yakovlev 
(1980) as a 2D integral which depends parametrically on 
the relativistic parameter x. Lampe (1968) analyzed this 
function in the static screening approximation at cc <C 1 . 
The asymptotes of J were obtained by Lampe (1968) for 
x <C 1 at 2 / ^ 1 and y » 1, by Flowers & Itoh (1976) for 
2 / » 1 at any x, and by Urpin & Yakovlev (1980) for 2 / 1 

and X » 1. Timmes (1992) calculated J{y) numerically in 
the limit of x 3> 1. However, the unified expression of J{y) 
at T < Tp valid equally for relativistic and non-relativistic 
electrons has been absent. Note that the fit expression of 
J{y) obtained by Timmes (1992) (his Eq. (10)) is valid 
only at 2 / < 10 ^. 

We have calculated J numerically for a dense grid of x 
and y in the intervals 0.01 < x < 100 and 0.1 < y < 100 . 
The results are fitted by the expression: 


J = 


6 2 
5x^ bx"^ 


3(1-h 0.0741422)3 


X In 


2.810-0.810/12-h 22 \ tt' 


y 


6 ( 13 . 91 -^ 22 )"^ 


(A3) 


which reproduces also all the asymptotic limits mentioned 
above. The mean error of the fits is 3.7%, and the maxi¬ 
mum error of 11 % takes place at x = 1 and y = 0 . 1 . 


A.2. Coulomb logarithms 

The tables of the Coulomb logarithms calcul ated from 
Eqs. & and (15) as described in Sect. ^.3.3| are htted 
by 


A.l. Electron-electron collisions 

The effective collision frequency Vee of non-relativistic de¬ 
generate electrons (x 1 , T < Tp) was analyzed by 
Lampe (1968) using the formalism of the dynamic screen¬ 
ing of the electron-electron interaction. The expression of 
t'ee for the relativistic degenerate electrons at T <C Tpe 
was obtained by Flowers & Itoh (1976). Here, Tpe is the 
electron plasma temperature determined by the electron 
plasma frequency Wpe, 

Tpe = hujpe/kB, UJpe = \/ He/ (Al) 

ml = me VT -b x2. The degeneracy temperature Tp in 
the NS envelopes is typically higher than Tpe. Urpin & 
Yakovlev (1980) extended the results of Flowers & Itoh 


I -b Cybe 


L = - (1 -b C5 &) In 


“ ^ /3^(1 - cef)) -b ^nZafd, 
where b = be + bi and 


- C4 


k = 


\ 2/3 


2TrZ. 


I.5ci 


C 3 VZ -b C 2/32 

Vtz 



For Z < 3, the fit parameters are given by: 
Cl = -0.349 -b 1.766 Z - 0.413 
C 2 = 3.499 - 2.355 Z -b 0.871 
C 3 = 11.92- 6.77Z-b 1.04 2-2, 

C 4 = 0.5733 - 0.1116 Z -b 0.0159 ^ 2 , 

C 5 = 0.7060 - 0.9896 Z -b 0.3517 ^ 2 , 


(A4) 
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C 6 = 4.779 - 2.676 Z + 0.361 
C 7 = 2.457 + 2.523 Z - 1.014 

For Z > 3, we obtain ci = 1.232, C 2 = 4.273, C 3 = 
0.9742, C 4 = 0.3816 C 5 = 0.9025, cg = 0.0, cy = 0.9. 

If Z = 1, the rms fit error over all our data set is 
about i5av ~ 3.2%, and the maximum error of ^max ~ 
7.2% takes place at p « 10^ g cm“^ and F = 3. If Z = 2 
we have Sa,v ~ 3.7%, iJmax ~ 7.8% at p « 10® g cm“^ 
and F = 1. For Z = 3, we obtain (5av ~ 2.4%, (Jmax ~ 
7.0% at log(p[g cm“^]) « 2.5 and F = 1. For higher Z, 
the rms error remains about 1.5-1.7%, and the maximum 
error mainly decreases. For instance, we have i5max ~ 5.1% 
at logp ~ 2.25 and F = 4 for Z = 6 ; <5max ~ 3.7% at 
logp ~ 5.75 and F = 2 for Z = 12; and (5max ~ 3.4% at 
logp « 5.5 and F = 8 for Z = 26. This fit accuracy is 
quite sufficient for studying the thermal structure of NSs. 

A.3. Relation between internal and effective temperatures 

In this section, we derive a fitting formula for Teff as a 
function of Tb, valid for 4.7 < log Teff < 6.5, 0.4 < 514 < 6 , 
where pi 4 is the surface gravity in units of 10 ^'^ cm s“®. 

Let us define Tbg = Tb/10® K, Tesg = Teff/lO® K, and 


■q = gl^AM/M, 


(A5) 


where M is the NS mass and AM is the accreted mass. 
According to GPE, 77 = Pb/(1-193 x 10®® Mbar), where 
Pb is the pressure at the bottom of the accreted envelope. 

For a purely iron (non-accreted) envelope, a very crude 
estimate (with an error ^ 30%) yields 


Teffe = T, = (7Tb9. (A6) 

Then C = Tbg — (T*/10^) is approximately an increase of 
the temperature through the iron envelope (in 10® K). A 
more accurate fitting formula reads 

T.m,Fe = 9i4i7Cr-^^ + iCm^-^% (A7) 

The typical fit error of Tesg.Fe is about 2%, with maximum 
4.2%, over the T^s — g domain indicated above. 

For a fully accreted envelope, we obtain 


if;Va = '?14(18.1Tb9)'-4®, 


(A 8 ) 


which is valid at not too high temperature (Tb < 10® K). 

Finally, for the partially accreted envelopes at any tem¬ 
peratures within the indicated range, we have 


rp4 

^em 


®i^e1l6,Fe 


rp4 

^ eff6,a 


d “h 1 


(A9) 


where 


a = [1.2 -b (5.3 X 10-®/^)°'^®] 


The typical fit error of Eq. (^) for Teffg is about 3%, with 
maximum 5.2%, for all possible values of q and any values 
of g and Teff within the indicated range. 

The dependence (A7) is recovered not only at suffi¬ 
ciently low accreted mass (p ^ 0 ), but also at sufficiently 
high Tb. The latter result reflects the fact (first demon¬ 
strated by GPE) that at high Teg the thermal insulation 
is mostly produced by the conductive opacities in the deep 
and hot layers of the envelope, in which light elements (H, 
He) burn into heavier ones. On the other hand, even at 
very low accreted mass (logAM/M ^ —16), the approx¬ 
imation of fully accreted crust is good enough at suffi¬ 
ciently low temperature because in this case the thermal 
insulation is provided mainly by the low-density accreted 
surface layers. 
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